/// @file qdist.cc
/// @brief Main program for quartet distance computation

// *************************************************************************

// Copyright© 2010 Texas A&M University, 
//                 College Station, Texas
//                 Contact: R Crosby (rwc@cs.tamu.edu)
//
// Copyright© 2003, 2004 by BiRC -- Bioinformatics Research Center
//                 University of Aarhus, Denmark
//                 Contact: Thomas Mailund <mailund@birc.dk>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

// *************************************************************************

#include <fstream>
#include <iomanip>
#include <math.h>
#include <sstream>
#include <string.h>
#include <time.h>

using namespace std;

#include "config.hh"				  // Autogenerated include needs to preceed other user includes

#include "BuildFlags.hh"
#include "dag.hh"
#include "options.hh"
#include "parser.hh"
#include "S.hh"

// *********************************************************************

static unsigned ParseCommandLine ( int                     argc,
				   const char *            argv[],
				   vector<const char*> &   tree_files );

static unsigned LoadDAG          ( vector<const char*> &   tree_files,
			 	   DAG &                  dag,
			 	   SMap &                 smap );

static void     PrintDAG         ( DAG &  dag,
			 	   SMap & smap );

static void     PrintStats       ( DAG &  dag );

// *********************************************************************

#ifdef HAVE_LIBPOPT

#include <popt.h>				  // The following include is here since is must follow config.hh

    static struct poptOption info_options[] = {
	{ "version", '\0', POPT_ARG_NONE, &options::version, 0, "Print version of qdist.",    0 },
	{ "contact", '\0', POPT_ARG_NONE, &options::contact, 0, "Print contact information.", 0 },
	{ 0 } // sentinel
    };

    static struct poptOption main_options[] = {
	{ "bare",       'b',  POPT_ARG_NONE,          &options::bare,       0, 
	  "Bare matrix output.",                   0 },
	{ "verbose",    'v',  POPT_ARG_NONE,          &options::verbose,    0, 
	  "Toggle verbose output.",                   0 },
	{ "ansicolor",  'a',  POPT_ARG_NONE,          &options::ansiColor,  0, 
	  "Turn off memorization of polynomial products"
	  " (exchanging runtime for memory usage).",  0 },
	{ "tree-limit", 'l',  POPT_ARG_INT,           &options::tLimit,     0, 
	  "Read no more than the specified number of trees "
	  "from the input files.",                    0 },
	{ "dissimilar", 'd',  POPT_ARG_NONE,          &options::dissimilar, 0,
	  "Process trees not all containing the same set of taxa", 0},
	{  0,           '\0', POPT_ARG_INCLUDE_TABLE, &info_options,        0, 
	   "Information options",                      0 },
	POPT_AUTOHELP
	{ 0 } // sentinel
    };

#endif

namespace {

    PerfVar pf("Mainline");			  // Performance trace block

}

// *********************************************************************
inline large_int_t
total_q_number(unsigned long n)
{
    // 2*binom(n,4)
    large_int_t res = n;
    return (res*(res-1)*(res-2)*(res-3)) / (4*3);
}

// *********************************************************************
int
main(int argc, const char *argv[])
{

    PRF(pf, "InitPgm");
    vector<const char*> tree_files;

    //
    // Parse the command line
    // 

    if ( unsigned rc = ParseCommandLine( argc, argv, tree_files ) != 0 ) 
	return rc;

    //
    // Process through all the input files loading the dag
    //

    SMap    smap;
    DAG     dag(smap);

    if (options::dissimilar)
	dag.Dissimilar(true);


    PRF(pf, "LoadDAG");
    if ( unsigned rc = LoadDAG ( tree_files, dag, smap ) != 0 )
	return rc;

    DB() cerr << dag << endl;

    //
    // Process the dag
    //

    PRF(pf, "CountDAG");

    dag.MakeHashTable ( smap.Size() );		  // Init the hash table
    dag.ClearMaps     ();			  // Clear dag loading maps
    dag.count         ();			  // Do the counting....

    //
    // Print the results
    //

    PRF ( pf, "Print" );
    PrintDAG ( dag, smap );

    //
    // Output stats
    //

    PrintStats ( dag );

    PRF(pf, "Final");
      
    return 0;
}

/// *********************************************************************
/// 
/// Process the command line arguments
///
/// \param argc Number of command line arguments
/// \param argv Array from main
/// \param tree_files vector of file names found
/// \return 0 if OK, value to return to os otherwise
///

static unsigned 
ParseCommandLine ( int             argc,
		   const char *    argv[],
		   vector<const char*> & tree_files )
{

#if HAVE_LIBPOPT

    poptContext ctxt = poptGetContext(0, argc, argv, main_options, 0);
    poptSetOtherOptionHelp(ctxt, "tree1 tree2 ..");

    int opt = poptGetNextOpt(ctxt);
    if (opt < -1) {
        cerr << poptBadOption(ctxt, POPT_BADOPTION_NOALIAS)
             << ':' << poptStrerror(opt) << endl;
        return 2;
    }

    if (options::version) {
        cout << PACKAGE_STRING << endl;
        return 0;
    }

    if (options::contact) {
        cout << PACKAGE_STRING << "\n\n";
        cout << "For questions or comments on qdist, contact:\n"
             << "\tRalph Crosby <rwc@cs.tamu.edu>\n";
        return 0;
    }

    while ( char *fname = (char*)poptGetArg(ctxt) )
        tree_files.push_back(fname);
    
#else
    for (int i = 0; i < argc; ++i)
	tree_files.push_back(argv[i]);
#endif

    return 0;
}

/// *********************************************************************
/// 
/// Load the trees from the input files
///
/// \param tree_files  Vector of the file names to load
/// \param dag         Reference to the DAG to load
/// \param smap        Reference to the S mapping block to use
/// \return 0 if no errors, value to return to OS otherwise

static unsigned 
LoadDAG ( vector<const char*> & tree_files,
	  DAG &                 dag, 
	  SMap &                smap )
{

    vector<const char*>::const_iterator f_itr;

    for ( f_itr = tree_files.begin(); 
          f_itr != tree_files.end() && 
	      (options::tLimit == 0 || dag.Roots().size() < options::tLimit); 
          ++f_itr) {

	try {
	    if (options::verbose) cerr << "Parsing " << *f_itr;
	    Parser::parse_tree( *f_itr, 
                                dag, 
			        smap );
	    if (options::verbose) cerr << endl;
	}
	catch (Parser::FileOpenError ferr) {
	    cerr << "Could not open file "
		 << ferr.fname() << endl;
	    return 1;
	}
	catch (Parser::ParseError perr) {
            cerr << endl 
		 << "Parse error in line " << perr.line_no()
                 << " in file " << perr.fname()
                 << endl;
            return 2;
	}
	catch (Parser::DegreeError derr) {
	    cerr << endl 
		 << "The tree in file " << *f_itr
		 << " is not binary.\n"
		 << "Currently only binary trees are accepted."
		 << endl;
	    return 3;
	}

    }

    if (dag.Roots().size() < 2) {
        cerr << "ERROR: at least two trees expected.\n";
#if HAVE_LIBPOPT
        cerr << "Use --help for details.\n";
#endif
	return 2;
    }

    if ( dag.Dissimilar() )
	cerr << "Trees containing disimilar taxa found" << endl;

    return 0;

}

/// *********************************************************************
/// 
/// Print the dag matrix
///
/// \param dag Reference to the dag to print

static void 
PrintDAG ( DAG & dag,
	   SMap & smap )
{

    // max_lineno is the printing width of the line number field

    stringstream        max_out;
    max_out << dag.Roots().size();
    const unsigned int  max_lineno      = string(max_out.str()).length();

    for ( unsigned int i = 0; i < dag.Roots().size(); ++i ) {

	if ( ! options::bare )
	    cout << setw ( dag.LabelWidth() )
		 << string(dag.Fname(i)).substr(0,dag.LabelWidth())
		 << "(" << setw(max_lineno) <<  dag.LineNo(i) << ")"
		 << " : ";

	// What results contains is the upper right triangular portion of the
	// actual full matrix with the exception that there may actually be
	// more entries in a row than are necessary. 
	// 
	// Only the upper right entries are actually used but the whole
	// matrix is printed...
	// 
	// The offset applied to the k value is to handle the extranious 
	// entries at the beginning of a row 

	for ( unsigned int j = 0, k = i; 
	      j < i; 
	      j++, k-- ) {

	    unsigned rowoff = dag[j].size() - ( dag.Roots().size() - 1 - j);

	    large_int_t prt = ( 
		                total_q_number(dag.Taxa(i)) + 
                                total_q_number(dag.Taxa(j)) - 
                                2 * dag[j][k - 1 + rowoff] 
 		              );

	    cout << prt / 4 << ((prt/2 & 1) ? ".5 " : " ");

	}

	cout << 0 << ' ';

	// The starting point for j is computed because there may be extranious 
	// entries at the start of a row

	for( unsigned int j = dag[i].size() - ( dag.Roots().size() - 1 - i),
	     l = i + 1; 
	     j < dag[i].size(); 
	     j++, l++) {

	    large_int_t prt = (
		                total_q_number(dag.Taxa(i)) + 
                                total_q_number(dag.Taxa(l)) - 
                                2 * dag[i][j]
                              );

	    cout << prt / 4 << ((prt/2 & 1) ? ".5 " : " ");

	}

	cout << endl;

    }

}

/// *********************************************************************
/// 
/// Output final statistics
///
/// \param dag Reference to the dag to process

static void 
PrintStats  ( DAG & dag  )
{

    if (options::verbose) {
	cerr << "DAG Statistics" << endl;
	cerr << "  Inner Nodecount called/computed:" 
	     << setw(11) << DagInfo::innerNodeCountCalled   << "/" 
	     << setw(11) << DagInfo::innerNodeCountComputed << endl;
	cerr << "  Inner Count called/computed    :" 
	     << setw(11) << DagInfo::innerCountCalled       << "/" 
	     << setw(11) << DagInfo::innerCountComputed     << endl;
	cerr << "  Target Traversals              :" 
	     << setw(11) << DagInfo::targetTraversals       << endl;
	cerr << "  Node/Vector memory used        :" 
	     << setw(11) << DagInfo::nodeMemory             << "/" 
	     << setw(11) << DagInfo::maxVectorMemory        << endl;
	cerr << "  Vector capacity                :" 
	     << setw(11) << DagInfo::maxVectorCapacity      << endl;
	cerr << "  Max vectors allocated          :" 
	     << setw(11) << DagInfo::maxVectors             << endl;
	cerr << "  Root/Inner/Leaf nodes          :"
	     << setw(11) << DagInfo::rootNodes              << '/'
	     << setw(11) << DagInfo::innerNodes             << '/'
	     << setw(11) << DagInfo::leafNodes              << endl;

    }

    DB(){
	// Get the vector distribution and print

	cerr << "DAG Vector Memory: " << endl;

	vector<unsigned> v = dag.VSizeDist();
	unsigned long mem = 0;

	for (unsigned i = 0; i < v.size(); i++) {
	    if (v[i]) {
		cerr << setw(5) << i << "," << setw(5) << v[i] << "," << setw(9) << v[i] * i << endl;
		mem += v[i] * i;
	    }
	}

	cerr << "Total Vector Memory: " << mem * sizeof(large_int_t)<< endl;
    }

}
